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Abstract 

Regular arrays of carbon nanotubes (CNTs) are frequently used in studies on field emission. However, non-uniformities 
are always present like dispersions in height, radius, and position. In this report, we describe the effect of these 
non-uniformities in the overall emission current by simulation. We show that non-uniform arrays can be modeled 
as a perfect array multiplied by a factor that is a function of the CNTs spacing. 
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Background 

Carbon nanotube (CNT) arrays for field emission (FE) 
applications have been extensively studied experimentally 
and theoretically [1-5]. Various improvements to fabricate 
well-aligned CNT arrays have been achieved, but non- 
uniformities are always present. To build precise arrays 
is expensive and difficult in extending to large areas. 
Simulation of CNT arrays is cost effective; however, simu- 
lation of these structures including non-uniformity is rare 
in the literature. To model non-uniformities in FE, it is 
necessary to understand their effects on the emission 
current. The simulation of FE in large domains is notori- 
ously difficult especially in three dimensions, which is neces- 
sary in this analysis. The difficulties include long simulation 
times, large computer memory requirements, and computa- 
tional instability. The first analysis of this kind is the recent 
work of Shimoi and Tanaka [6]. They managed to perform 
three-dimensional (3D) simulations based on boundary ele- 
ments that avoided meshing the volume of the 3D domain. 
They simulated carbon nanofibers with random position 
and height to match the emission pattern that they obtained 
experimentally. In this work, we perform simulations of 
non-uniform CNTs with dispersions in height, radius, 
and position in limited ranges and with small CNT aspect 
ratios aiming to correlate the current from non-uniform 
arrays with the current expected from perfect arrays. 
We restrict our analysis to a hemisphere-on-a-post model 



[4,6-8], in which the CNTs are regarded as perfect con- 
ductors, with a smooth surface and oriented normal to 
the substrate. In this report, we shall refer to these ide- 
alized tubes as CNTs. 

Methods 

The CNTs are positioned in a 3 x 3 square array, as shown 
in Figure 1. We shall explain hereafter that a 3 x 3 square 
array is an efficient way to perform the simulations. The 
ith CNT height H b radius R it and coordinates (X it Yd are 
stochastic variables with expected values (or averages), re- 
spectively, equal to h = 10 a.u., r = 1 a.u., and (x b yi) being 
the center of the /th unit cell in the array. Thus, the de- 
fault aspect ratio is 10, which is quite small. However, 
larger aspect ratios cause simulation difficulties that will 
be discussed later. Despite this limitation, we think that 
the results provide a meaningful insight on the behavior 
of the current. We simulated aspect ratios up to 100 in 
graphenes randomizing only the positions. The results 
vary at most 25%, tending to increase slowly in a loga- 
rithmic pace as a function of aspect ratio. A complete 
analysis of graphene sheets will be presented in a forth- 
coming paper. The stochastic variables in our study will 
be limited to the following ranges: 
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Figure 1 Hemisphere-on-a-post model for a 3x3 non-uniform array domain. In (a), (b), and (c), respectively, height, radius, and position 
are separately randomized. In (d), all three parameters are randomized at the same time. The red regions indicate strong electric field. 
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where s is the array spacing; a h) a n and a p can be 
interpreted as the range in percentage of the expected 
value. For instance, a h = l implies that the height of the 
CNT can vary 100%, from 0.5 h to 1.5 h. The choice for 
these dispersion ranges was based on microscopic obser- 
vations [6,9,10]. If a = 0, the corresponding stochastic vari- 
able is constant. Equation (3) states that the displacement 
range of the CNTs can vary from no displacement (a p = 0) 
to displacements as large as half the length of the unit cell 
(a p = l). We analyze the emission current as a function of 
s from near close packed (s > 0.25 h) to s = 10 h (approxi- 
mately isolated tubes). The field enhancement and the 
screening effects are illustrated in Figure 1. In Figure la, 
only the heights are randomized. The taller the tube, the 
larger the field strength at the tip, represented in shades of 
red; shorter tubes are shielded. In Figure lb, only the radii 
are randomized. The screening effect is approximately the 
same for all tubes, but the field enhancement is larger at 
the thinner ones. In Figure lc, only the positions are ran- 
domized. In this case, some tubes are more screened than 
others depending on how they clump up, notice however, 
that the field strength at the tips are more homogeneous 



compared to Figure la,b. Indeed, the overall current is less 
affected by randomized positions than heights or radii for 
the separation shown in this figure. In Figure Id, all vari- 
ables are randomized at the same time. The CNTs are not 
allowed to overlap. 

The simulations are performed using software COMSOL® 
v.4.2a, which is based on the finite elements method. The 
CNT array, as shown in Figure 1, is regarded as purely elec- 
trostatic system. A macroscopic vertical electric field of 
10 GV/m is applied on the domain. The side boundaries 
have symmetry boundary condition, which states that 
there is no electric field perpendicular to these boundaries 
(E.n = 0) making them act as mirrors. These conditions 
determine the norm of the electric field in the domain. 

The local current density, is evaluated using Fowler 
Nordheim equation [11,12]: 

. AE 2 ( B0 3 /2\ 

where A = 1.56 x 10" 6 A eV V" 2 , B = 6.83 x 10 9 eV" 3/2 V/m, 

0 is the work function (in eV), and E is the local electric 
field (in V/m) at the surface of the CNTs. We use a work 
function of 5 eV for the CNTs. Equation (4) is integrated 
over the CNT's surfaces to obtain the overall current, 
which is normalized by the current from a perfect array 

1 pa* Figure 2 shows I PA and the overall current density, 
Jpa, defined as the total current divided by the area of the 
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Figure 2 Field emission current l PA and current density J PA of a 
perfect array. The lattice spacing s is expressed in units of the CNT 
height h. The aspect ratio of the CNTs is 10 in this figure. 



array. The peak in J PA a.ts = 2h indicates the ideal spacing 
for FE applications [13,14]. Note that J PA is relatively small 
for s < h, so we shall focus most of our analyses to the 
region where s > h. The currents and current densities 
shown in Figure 2 for the perfect uniform lattice and 
uniform CNTs will be used to normalize the currents 
for the non-uniform structures. 

Each simulation run, identified with the number of the 
run, /<, has a particular set of randomized parameters 
that yield the normalized current, 4- The I k values from 
a 3 x 3 domain present large variations, but after aver- 
aging 25 simulation runs, we obtain a smoother behavior, 
which is the expected values of the stochastic 4 The error 
in 4 decreases by a factor of lHk. In FE experiments, the 
observed current is the average over a large number of 
CNTs. We did 25 simulation runs of 3 x 3 CNTs, which is 
physically similar to simulate 225 CNTs in one run. How- 
ever, the latter calculation is impossible due to memory 
and numerical instability. Even a 3 x 3 array takes a rather 
long time to simulate, and some of our results were not 
reliable at large spacing. We simulated arrays with lxl, 
2x2, 3x3, and 4x4 randomized CNTs. The average 
current depends on the size of the domain, but the con- 
vergence is fast. The normalized currents as a function of 
the spacing for 3 x 3 and 4x4 arrays are exactly the same 
within the error. Hence, a 3 x 3 domain is already large 
enough to represent a random field of CNTs. 

Results and discussion 

Figure 3 shows the result when only the positions of the 
CNTs are randomized (a p =1, a r = a h = 0). The normal- 
ized average I p = <4 > is shown in full circles. The gray 
line at I p = 1 is drawn to guide the eye. The sine-like be- 
havior of I p is a consequence of the step shape of I PA 
(see Figure 2), which increases fast at small s and satu- 
rates for s — > oo. The random positioning causes some 
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Figure 3 Randomization in the (x,y) coordinates of the CNTs in 
the array. The gray opened circles are the normalized current l k 
from an individual simulation run. The full circles are the average 
over 25 runs (/ p ). The inset shows s> h superposed to an interpolating 
function that provides a numerical value for \ p . 



CNTs to lump, while others form a sparser configur- 
ation. At small 5, the field enhancement of the slightly 
isolated CNTs dominates the lumping of CNTs elsewhere, 
thus I p > 1. On the other hand, for large s, the CNTs 
are practically isolated, and their field enhancement of 
the CNTs is almost at a threshold value. In this case, 
the current from isolated CNTs is almost constant, while 
the screening effect of the lumped regions significantly 
reduces the current, so I p < 1. For s — > oo, the emitters 
are isolated, and it is unlikely that two or more emitters 
will become close enough to screen each other after 
random displacements; therefore, I p tends to unity. At 
s = h, field enhancement and screening on the random- 
ized tubes compensate exactly and I p = l. At this point, 
misplaced CNTs do not affect the overall current expected 
from a perfect array. The inset in the figure shows the 
region for s > 1, which is the important region for FE ap- 
plications as mentioned. We fitted this region with the 
simplest interpolating function to provide a numerical 
value for I p . The fitting curve is shown in the inset. 

Figures 4 and 5 show the normalized currents I r and I h 
for a r = 1 and a h = l, respectively. Like in Figure 3, the 
horizontal axes in these figures are logarithmic. At small 
5, I n and 4 are sensitive to the randomization as can be 
seen. In this region, fluctuations in height and radius 
largely decrease the electrostatic shielding as compared 
to the uniform CNTs, thus the normalized current be- 
comes very high. It should be remembered that, although 
the normalized I r and I h are high for small 5, the absolute 
current is actually very small, as can be seen in Figure 2. 
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Figure 4 Normalized current from randomized radii of the CNTs. 
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The insets show the curves for s > h. The interpolating 
functions used in Figures 3, 4, and 5 for s > h are 

I p (s >h) = 1.02-1.9(s/h + 1.01)" L845 ln(5//z), (5) 
I r (s >h) = 1.26 + 17.7(s/h + 0.833)" 4 - 884 , (6) 

I h {s >h) = 1.41 + 579{s/h + l)" 8766 + 0.04771n(s//z). 

(?) 

Equations (5) to (7) have no physical meaning; they 
are mere interpolating functions only to provide numer- 
ical values between the simulated points. These inter- 
polating functions were chosen for representing the 
shape of the curves by taking the logarithmic scale of 
the #-axis into account. 
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Figure 5 Normalized current from randomized heights of the CNTs. 
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Next, we analyze the effect of randomizing two parame- 
ters simultaneously. It is not trivial to evaluate, for example, 
I pr knowing the values of I p and I r . The difficulties are the 
non-linearity of Eq. (4) and the complicated local electric 
field E that appears in it. This field is a function of X b Y b R t 
and H t and does not have an analytic solution. Therefore, 
for this analysis, we need to vary two parameters simul- 
taneously. Just as for I p , I r or I h , the simulations are av- 
eraged over 25 runs. The results are shown in Figure 6. 
In this figure, the expected values of the normalized 
current are specified with two sub-indices that indicate 
the parameters that are varying. Figure 6 also shows the 
expected normalized current I prhi when varying the 
three parameters: position (x,y), radius, and height at 
the same time. Interestingly, I prh is below the curves for 
Ihr and I p h in some regions. This means that randomiz- 
ing two parameters affects the average current more 
than varying three parameters in these regions. The 
curves are always greater than unity, typically between 1 
and 4 for s > h. This is a consequence of randomization: 
some CNTs are less electrostatically screened causing them 
to surpass the emission of a perfect array. Furthermore, 
most CNTs are screened, as can be seen in Figure Id; so, 
only few CNTs are accounting for the total current [6]. 
Then, by increasing the external electric field, these few 
CNTs will become overloaded before most CNTs can start 
contributing to the current. Consequently, the maximum 
current density of non-uniform arrays is limited by the 
current that these few CNTs can support. We define 
^high as the highest CNT normalized current in the 3x3 
array averaged over 100 runs. 7 h i g h comprehends 1/9 or 
11.1% of the most emissive CNTs. Figure 7 shows 7 h i g h 
as a function of s for s > h and its standard deviation, 
0^high> shown in the figure as error bars. The a-/ h i g h can 
be used to determine what part of the CNTs is expected 
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Figure 6 Normalized emission randomizing variables two at a 

time and all three variables simultaneously. 
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Figure 7 Highest normalized emission / high and the standard 
deviation o7 h j g h as a function of the spacing. The a/ hi gh is shown 
as half error bars. These parameter can be used to estimate the 
fraction of CNTs that will burn out at certain current given the 
degree of non-uniformity. 



to burn in the non-uniform array given their tolerance, 
as we shall indicate below. 
The interpolating functions for the curves of Figure 6 are 



where I pr is given in Eq. (9). Finally, if all a parameters 
vary, we have 

l(a p ,a r ,ah) =l(a p ,a r ,6) + a 2 h (l phr -l(a p ,a n 0)), 

(14) 

where I phr is given in Eq. (11). 

From the data shown in Figure 7, we derive the follow- 
ing interpolating functions 

4igh (s > h) = 1 + 1.12V/* + 3 - 34 </* + 19a ^ exp ("^^) ' 

(15) 

where, a prh = max(a p a r a h ) and 

/ s/h-1 

<r/ high (s >h) = 0.16a prh + 2.42^ + 16.1a 2 prh expl--^ 

(16) 

Equations (15) and (16) give an upper estimate of the 
maximum current carried by individual CNTs, as a func- 
tion of our randomization parameter a pr ^ 

The fraction of CNTs expected to burn out can be 
evaluated from a Gaussian distribution as: 



I ph (s >h) = 1.09 + 38.2(s//z + 0.6)" 6 ' 235 + 0.148ln(s//z), 

(8) 



Ipr(s >h)= 0.93 + 6.08(s//z + 0.72)" 319 + 0.091n(s//z), 

(9) 



I hr (s > h) = 1.75 + 2.68(s//z + 0.033) 



-4.21 



(10) 



I phr (s >h) = 1.31 + 0.5(s//2-0.23)" JM * + 0.28ln(s//z). 

(11) 

Equations (5) to (11) are valid for a=l; however, our 
simulation results (not shown here) indicate that a quad- 
ratic function fits intermediate values 0 < a < 1 reason- 
ably well. The following example gives a procedure to 
obtain the normalized current for any set {a pi a n ay), with 
normalized current I{a p ,a„ap). In the simplest example, 
if only a p varies, then 



7(^,0,0) =1 + 02(7,-1), 



(12) 



where I p is given by Eq. (5). In another example, in 
which a p and a r are varying, then 



I(a p ,a r ,0) =7(^,0,0) +^(7^-/(^,0,0)), (13) 



11.1% .„ / («'-7 high 



c7 h ighv^r 

±(l+erf( Ihiih ~ Ina A ) xll.1%, 



(17) 



where erf(z) is the error function, 7 max is the normalized 
burn out current (or tolerance). The factor 11.1% is be- 
cause Eqs. (15) and (16) account only for l/9th of the 
CNTs in the 3x3 array. 

Let us give an example: consider a non-uniform array 
with a p = 0.4, a r = 0.5, a h =0.8 observed microscopically 
and s = 2h yielding an average emission of 1 \xA. From 
Eqs. (14), (15), and (16), we calculate a normalized current 
of 7=1.28, which corresponds to the 1 \iA; 7 high = 4.94 
(3.86 \iA) and o7 high = 1.90 (1.48 \iA). Now, suppose 7 max is 
10 (7.81 |iA), then the fraction f of emitters that will burn 
out at 1 \iA is smaller than 0.04% according to Eq. (17). In 
this example, 7 max is constant: otherwise, the calculation 
of f will be more elaborate. If 7 max is a known function, 
then f must be integrated over 7 max for a refined estima- 
tive. However, we shall not deepen our analysis on f in 
this paper. 

Conclusions 

We simulated the behavior of the field emission current 
from non-uniform arrays of CNTs and obtained correc- 
tion factors to multiply the current from a perfect CNT 
array toward the currents of non-uniform arrays. These 
correction functions are valid if the allowed dispersion 
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in height and radius is kept inside the limits of 50% and 
150% of their average values and if the randomization of 
the CNT position is done inside the designated unit cell. 
The uneven screening effect in non-uniform arrays 
causes many CNTs to become idle emitters while few 
may become overloaded and burn out. To avoid this, 
uniformity is desired: however, non-uniformities are al- 
ways present in some degree, and our model describes 
how to treat them. This model can also be used in es- 
timating how many CNTs are expected to burn given 
their tolerance and the total current extracted from 
the array. 

We like to point out that in a previous work [15], we 
showed that the emission from 3D CNT arrays can be 
simulated in a two-dimensional (2D) rotationally symmet- 
ric system with proper boundary conditions. The currents 
from the 2D and 3D arrays are also related by a factor that 
is a function of the aspect ratio and spacing of the actual 
array. The combined correction factor from Eq. (14) and 
the procedure in [15] can considerably ease the modeling 
of FE from non-uniform CNT arrays, as they can be re- 
duced to perfectly uniform arrays, which may be treated 
in a 2D model. 
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